Thermal Contraction and Disordering of the Al(llO) Surface 
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Al(110) has been studied for temperatures up to 900 K via ensemble density-functional molecular 
dynamics. The strong anharmonicity displayed by this surface results in a negative coefficient of 
thermal expansion, where the first interlayer distance decreases with increasing temperature. Very 
shallow channels of oscillation for the second-layer atoms in the direction perpendicular to the 
surface support this anomalous contraction, and provide a novel mechanism for the formation of 
adatom- vacancy pairs, preliminary to the disordering and premelting transition. Such characteristic 
behavior originates in the free-electron-gas bonding at a loosely packed surface. 

PACS numbers: 71.15.Pd, 65.70.+y, 68.35. Ja, 68.45.Gd 
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Metal surfaces exhibit a remarkable behavior as a func- 
tion of the temperature. Thermodynamic stability is of- 
ten determined by a delicate balance between energetic 
and entropic effects, and can lead to a rich phenomenol- 
ogy for the phase diagrams of different systems. Unrecon- 
structed face-centered cubic (110) surfaces (e.g. Al, Cu, 
Ni) display a damped oscillatory pattern of interlayer re- 
laxations, starting with a large contraction between the 
first and the second layer |l]j2| . Such behavior originates 
in the response of surface atoms to under-coordination: 
moving towards the underlying layer, they increase their 
surrounding charge density while reducing the corruga- 
tion of the surface and the lateral tensile strain M . When 
the temperature is raised, this under-coordinated layer 
can start to disorder even before the melting temperature 
of the bulk is reached. While the suggestion that a sur- 
face could act as a nucleation stage for melting had long 
been made, experimental evidence of a reversible melt- 
ing transition limited to the outer surface layers came 
only recently H. For the case of Al(110), several experi- 
mental techniques (ion blocking and shadowing, electron 
or neutron diffraction, He scattering) have since shown a 
clear onset of disordering at temperatures between 770 K 
and 815 K whereas the bulk melting temperature is 
933 K. Computer simulations based on different models 
(effective-medium theory j^,f5J, embedded-atom method 
M , glue models [^) ) have then been applied to the study 
of several (110) surfaces (Pb, Al, Cu, Ni), and surface 
premelting was observed in all cases. 

However, many issues remain unresolved. Exten- 
sive low-energy electron diffraction (LEED) studies in 
Al(110) 0](c) show a negative thermal expansion coeffi- 
cient for the first interlayer distance, and a large positive 
one (twice the bulk value) for the second interlayer dis- 
tance. These findings are at variance with widely held 
general theoretical considerations |l(J, and with the re- 
sults of available computer simulations for Cu, Ni, and 
Al pl-fl which predict an expansion of the first interlayer 



distance with temperature. In addition, model calcula- 
tions fail to reproduce the zero-temperature multilayer 
relaxation pattern |?]||, predicting only the contraction 
of the first interlayer. On Al(110) the premelting transi- 
tion is preceded by an anomalous proliferation of adatoms 
on the surface 0, for which there is no reliable micro- 
scopic picture. Finally, the degree of anharmonicity and 
anisotropy of the different surface layers, as opposed to 
the bulk, is not known, due to the experimental difficulty 
in resolving different layers. 

Ab-initio molecular dynamics (MD) simulations of 
metal surfaces are very challenging, and only few and 
limited studies have been attempted. We use here an 
approach that we recently introduced jll],[l2| (ensemble 
density- functional theory (eDFT)), together with a tech- 
nical improvement for the Brillouin Zone (BZ) integra- 
tions (so called "cold smearing"), which is particularly 
suited to MD simulations. Applying this scheme to the 
case of Al(110), we provide the first theoretical confir- 
mation of a negative thermal expansion for this surface, 
and excellent agreement with the experiments for the 
temperature-dependent multilayer relaxations. More- 
over, we present a novel, and in retrospect simple, picture 
of the microscopic mechanisms that lead to this anoma- 
lous thermal contraction and to the surface disordering 
associated with premelting. 

In first-principles calculations for metals it is custom- 
to introduce a fictitious electronic temperature a 
to broaden the density of states and to smooth 
the discontinuities at the Fermi energy fi, greatly improv- 
ing the sampling accuracy of a given set of k-points. It 
is very convenient to choose a broadening that has zero 
first- and second-moments, so that the resulting elec- 
tronic free energy does not have any quadratic depen- 
dence on the broadening temperature, and neither do its 
derivatives with respect to any external parameter jl2 15 
(e.g. the Hellmann-Feynman forces, or the stress ten- 
sor). In the existing schemes this is achieved at the price 
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of allowing for negative orbital occupancies JHJ , so that 
problems can arise in self-consistent calculations where 
the total electronic density may become negative. Here, 
we present a broadening scheme leading to an occupation 
function that is positive definite. Occupation broadening 
convolutes the density of states with a broadening of the 
5 function Jl3|,[l4j ; the cold-smearing broadening is 
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Spin-degeneracy is assumed here, and x — U^l, The 
"generalized entropic functional" S = Si that can 
be derived |l3| ] and the occupation numbers /,; — 
$( x ) dx can all be expressed in terms of pseudoen- 
ergies ei |l(| (xi — !± -^ ± )', in particular 
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No practical difficulty to the self-consistent calculations is 
caused by the fact that some spin-degenerate occupancies 
can still exceed 2; this was also the case for the choices 
set forth in Ref. |l5). 

The calculations use the local-density approximation 
(LDA) and norm-conserving pseudopotentials, with a 
plane- wave basis cutoff of 11 Ry. The bulk proper- 
ties of Al are well represented: the lattice parameter is 
3.96 (4.02) A, the elastic constants C u =117 (114) GPa, 
Ci2=66 (62) GPa, and C 44 =39 (32) GPa (experimental 
results at K Jl7| are in parenthesis). The simulation 
cell is a 3 x 3 8-layer Al(110) slab, containing 72 atoms 
separated by 8.5 A of vacuum, k-point sampling is per- 
formed with the 4, i, i Baldereschi point, using o~=0.5 
eV of cold smearing. 

The zero-temperature structural properties are sum- 
marized in Table Q: good and consistent agreement with 
the experimental results is registered. The 8-layer calcu- 
lation has been performed using the same finite cell and 
sampling of the MD simulations; this introduces some 
small finite-size errors, that can be evaluated exactly at 
K comparing them with a fully converged calculation 
(a 1 x 1 15-layer slab with 12 x 12 x 2 k-point sampling). 

Constant-temperature MD simulations have been per- 
formed using a Gaussian thermostat and a leapfrog ve- 
locity Verlet algorithm to integrate the ionic equations 
of motion 18 1, using a timestep of 8 fs. A set tolerance 
for each timestep of 5 meV/cell in the spread of the to- 
tal energies over the last 5 electronic iterations resulted 
in a negligible drift of the constant of motion (less than 
1.5 meV/atom/ps for a microcanonical run). The lattice 
parameter parallel to the surface was fixed applying the 
experimental thermal expansion coefficient for the bulk 
to the LDA equilibrium lattice parameter. We followed 
five runs, at increasing temperatures of 400, 600, 700, 800 
and 900 K, for 5, 10, 6, 6, and 6 ps respectively. 
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FIG. 1. Layer-resolved mean square displacements at 
400 K; left panel shows the [001] components (henceforth la- 
beled x), the right panel the [110] ones (z). The third- and 
fourth-layer data are given by the unlabeled solid and dashed 
thin lines. The shaded area corresponds to the bulk experi- 
mental values of Ref. E5] . The horizontal axis shows the time 
over which averaging was performed. 

We present in Fig. [l] our results for the mean square 
displacements (MSDs) in the different layers, from the 
surface to the interior of the crystal, during a 5 ps run at 
400 K. The horizontal scale shows the decremental time: 
the plot starts with the averages over the full run, then 
proceeds by discarding a progressively longer initial seg- 
ment. This approach highlights the initial thermalization 
time (negligible in this case), the flatness of the plateau 
for the converged time average, and provides an estimate 
for the statistical errors. The left panel of Fig. |l| shows 
the [001] component for the MSDs (we label it x); this 
component is parallel to the surface and perpendicular to 
the [110] rows that characterize the (110) surface. The 
right panel shows the [110] z component, perpendicu- 
lar to the surface. The time averages are well converged, 
with the third- and fourth-layer results very close to each 
other (giving us confidence on the absence of finite-size 
effects), and close to the experimental bulk values |]l9| . 

Two results stand out from the simulation. First, the 
MSDs in the x direction are twice as large for the surface 
atoms than for those in all the other inner layers. While 
it can be expected that the undercoordinated atoms on 
the surface should be more loosely bound, the large dif- 
ference with the averages for the lower layers is notable. 
Second, the MSDs in the z direction (i.e. perpendicular 
to the surface) are much larger in the second layer than 
in the first layer. This is a distinctive feature of this 
crystallographic orientation that was first encountered in 
embedded-atom simulations of Ni(110) and Cu(110) [||. 
In Al the effect is more striking due to its free-electron- 
gas behavior. A simple rationalization can be offered: 
since the (110) surface is very open, atoms in the second 
layer have natural channels of oscillation perpendicular 
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FIG. 2. Layer-resolved mean square displacements at 
700 K. The x,y ones parallel to the surface are in solid black 
and dashed grey (x is perpendicular and y is parallel to the 
[110] rows); the z ones, perpendicular to the surface are in 
dashed black. "Bulk" refers to the two inner layers. 

to the surface and directed towards the vacuum. The 
charge density on the top of the second-layer atoms is still 
quite homogeneous, and the bonds are easily stretched, 
leaving thus the freedom for the atoms to move back and 
forth along these channels. On the other hand, atoms in 
the first (surface) layer see the vacuum acting as a hard 
wall, limiting their mobility outwards; their largest oscil- 
lations are thus parallel to the surface and perpendicular 
to the [110] rows. 

The anisotropic behavior of the surface dynamics can 
be gauged by looking at Fig. |^, where the MSDs at 700 K 
are plotted in all three crystallographic orientations as a 
function of the layer depth. Moving from the bulk to the 
surface, one can observe that the third layer still behaves 
in a bulk-like fashion: the MSDs are isotropic, and they 
are only slightly larger than those in the two layers be- 
low. The anisotropy becomes very distinct in the second 
layer, with its characteristic large MSDs perpendicular 
to the surface, and persists in the first layer, for which 
the 'easy' channels are parallel to the surface and across 
the close-packed [110] rows. Some of the components for 
the MSDs in the first and second layers can be up to 
2-3 times their bulk counterparts. These enhancements 
near the surface are due to the lower coordination; in 
addition, a higher degree of anharmonicity makes these 
surface MSDs increase much more rapidly with temper- 
ature than the bulk ones. This becomes apparent from 
the plot in Fig. || of the MSDs as a function of the tem- 
perature. The innermost layers show isotropic MSDs, 
with some deviation from the linear regime only above 
700 K (in the harmonic regime the MSDs dependence 
on the temperature is exactly linear). The outer lay- 
ers, on the contrary, are strongly anharmonic. The very 
large increases in the vibrational amplitudes along the 
'easy' channels are precursors to the creation of adatoms 
and vacancies on the surface, that lead to the disordering 
and premelting of the surface. In fact, we observe that 
with increasing temperature atoms in the second-layer 
start making increasingly large slow excursions towards 
the surface. One of these events is shown in Fig. |]; the 
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FIG. 3. Layer-resolved mean square displacements as a 
function of temperature. The shaded area corresponds to the 
bulk experimental values of Ref. |J19| . 
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FIG. 4. Projected view of the MD trajectories in the (110) 
plane, for the top 4 layers of the slab. The large excursion to- 
wards the surface of a second-layer atom has been highlighted. 

highlighted atom temporarily pops out from the sur- 
face. In another event the second-layer atom remained 
outside the surface, creating an adatom-vacancy pair 
where the vacancy is initially in the second layer (this 
void is quickly filled up by a surface atom) . In this second 
case the adatom diffused away via exchange diffusion. 

The microscopic dynamics provides a clear explanation 
of the behavior of this surface, that displays an increasing 
contraction of the first interlayer distance with tempera- 
ture, where a large expansion would have been expected. 
The contraction can be understood looking at the mo- 
tion of the second-layer atoms along these channels that 
are shallower towards the vacuum. With increasing tem- 
perature, the center of mass of the second layer moves 
outwards, since it is not hampered by nearest-neighbors 
directly on top (the first layer is staggered with respect 
to the second). The first layer is more limited in its ex- 
pansion, since the vacuum acts as a hard wall. The end 
result is that the average distance between the first and 
the second layer decreases with temperature. This de- 
crease is then offset by a larger thermal increase of the 
distance between the second and the third layer. The 
results of our simulations for the interlayer relaxation as 
a function of temperature (see Fig. |^) are in very good 
quantitative agreement with the LEED data (Ref. [0(c)). 
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FIG. 5. Interlayer relaxations as a function of temperature. 
Solid circles: eDFT MD simulations, with their statistical 
error bars; dashed lines: quadratic least-squares fit; crosses: 
experimental LEED values (Ref. [§(c)). The eDFT data 
have had their K finite-size corrections added rigidly at all 
temperatures. 

In conclusion, our calculations on Al(110) represent the 
first extensive first-principles molecular-dynamics simu- 
lations of the dynamics on a metal surface, presenting 
both an insightful picture of the microscopic dynamics 
and a remarkable agreement with the available experi- 
mental data. The microscopic dynamics of this surface is 
peculiar, and governed by the interplay between the free- 
electron-gas behavior of the bulk and the quasi-covalent 
bonding of the undercoordinated surface atoms. Two dis- 
tinct soft channels of oscillation have been identified. One 
channel is at the surface in the [001] direction, perpendic- 
ular to the close-packed surface grooves. The other, un- 
expected, is perpendicular to the surface but confined to 
the second-layer atoms. It is this channel that is responsi- 
ble for the observed anomalous contraction of the surface 
with temperature. Additionally, it provides a novel, fa- 
vored mechanism for the generation of adatom-vacancy 
pairs, whose proliferation is precursor to the disordering 
and premelting transition. 
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TABLE I. Theoretical and experimental values for the in- 
terlayer relaxations in Al(110); theoretical values are at K. 
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